function Master_Policy_EEGap_u_standard_Heterinc_co2_m_adj_pos_Wglow(subsample,inc,set_seed,matlab_seed,nb_draws)

% Created: 08/13/2018
% This code computes the optimal tax and welfare under three scenarios:
% 1. No adjustment for misperceptions: tax=damage
% 2. Adjustment taking into account the distribution of misperception
% 3. Adjustment taking into account the mean misperception

% This version superseeds v3. It uses the fixed effects from the mixed logit instead of the clogit.

% The results from 
% 
% Dependencies:
% Fox model v1 estimated for all income groups
% Master_Fox_mixed_discrete_FEdemoshrt_3param_grID_wtp_j_bs(44000,16,1,1,50, j, -)
												

% Parameters for the policy simulation

  elec_scale=1;  % Ensure that we keep heterogeneous electricity price
  p_elec=0;      % Ensure that we keep heterogeneous electricity price
% p_elec=0.11;   % Homogeneous electricity price
% kwh;
% elec;
% estar;
%  d_ext=0.02;   % Homogeneous carbon externality $/kEh
  d_ext_carbon_price=50; % Heterogeneous carbon externality $/ton CO2
  r=0.05;        % Market discount rate
  lifetime=18;   % Expected lifetime of the appliance
  MCPF=1;

  rho_5=1/(1+r);
  LLC_5_18=rho_5*(1-rho_5^18)/(1-rho_5);
  
  rho_2=1/(1+0.02);
  LLC_2_18=rho_2*(1-rho_2^18)/(1-rho_2);
  
  rho_12=1/(1+0.12);
  LLC_12_18=rho_12*(1-rho_12^18)/(1-rho_12);
  
  rebate_estar_denom_scale=0;

  Nsample=3500;

tic

%Create Diary
date_now = clock;
date_now = strcat(num2str(date_now(1)),'_',num2str(date_now(2)),'_', num2str(date_now(3)), num2str(date_now(4)), num2str(date_now(5)));

 %path_data = {'/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Data/Matlab_estimation/'};
 %path_results = {'/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_Simulation/'};

 %path_data = {'C:\Users\ecmyers\Dropbox\Appliance_EnergyPrice\EEgap_data_code_heter_SM\Data\'};
 %path_results = {'C:\Users\ecmyers\Dropbox\Appliance_EnergyPrice\EEgap_data_code_heter_SM\Results_PolicyAnalysis\'};

 path_data = {'C:\Users\MatlabGEM\Dropbox\eegap\EEgap_data_code_heter_SM\Data\Matlab_estimation\'};
 path_results = {'C:\Users\MatlabGEM\Dropbox\eegap\EEgap_data_code_heter_SM\Results_Simulation\'};


%Read Files
    %bchoice=load('H:\Research\sears\estar_data\refrigerators\bchoice_046_2008_2012_v11022017_struct_52171_4.csv'); 
    %tmp_file1=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\bchoice_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file1=strcat(path_data{1},'bchoice_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    bchoice=load(tmp_file1); 
    bchoice(:,1)=[];
    MT=size(bchoice,1);
    bchoice(Nsample+1:MT,:)=[];
    bchoice=bchoice';
    [bchoice_i,bchoice_j,bchoice_s]=find(bchoice);
    bchoice_ij=find(bchoice);
    [bchoice_m,bchoice_n] = size(bchoice);
%     clear bchoice;
    
    %tmp_file2=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file2=strcat(path_data{1},'choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    id=load(tmp_file2); 
    id(:,1:7)=[];
    id(Nsample+1:MT,:)=[];
    id=id';
    [id_i,id_j,id_s]=find(id);
    id_ij=find(id);
    [id_m,id_n] = size(id);
    
    
    MTvector=ones(MT,1);
    J=size(id,1);
   
    %tmp_file3=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\estar_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file3=strcat(path_data{1},'estar_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    estar=load(tmp_file3); 
    estar(:,1:7)=[];
    estar(Nsample+1:MT,:)=[];
    estar=estar';
    estar_num=estar(bchoice_ij);
    estar_denom=estar(id_ij);
    prob_estar_exact=mean(estar,1);
    %clear estar;
    
    %tmp_file4=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\rebate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file4=strcat(path_data{1},'rebate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    rebate=load(tmp_file4); 
    rebate(:,1:7)=[];
    rebate(Nsample+1:MT,:)=[];
    rebate=rebate/1000;
    rebate_J=repmat(rebate,1,J);
    rebate_J=rebate_J';
    rebate_denom=rebate_J(id_ij);
    rebate_estar_num=rebate_J(bchoice_ij).*estar_num;
    rebate_estar_denom=rebate_J(id_ij).*estar_denom;
    rebate_estar_denom=rebate_estar_denom_scale*rebate_estar_denom;
    clear rebate rebate_J;
    
    %tmp_file5=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\tax_rate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file5=strcat(path_data{1},'tax_rate_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tax=load(tmp_file5); 
    tax(:,1:7)=[];
    tax(Nsample+1:MT,:)=[];
    taxp=1+repmat(tax',J,1);
    
    %tmp_file6=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\price_mean_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file6=strcat(path_data{1},'price_mean_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    price_net=load(tmp_file6); 
    price_net(:,1:7)=[];
    price_net(Nsample+1:MT,:)=[];
    price_net=price_net';
    price=taxp.*price_net;
    price=price/1000;
    price_num=price(bchoice_ij);
    price_denom=price(id_ij);
    clear price_net taxp tax;
    
    %tmp_file7=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\kwh_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file7=strcat(path_data{1},'kwh_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    kwh=load(tmp_file7); 
    kwh(:,1:7)=[];
    kwh(Nsample+1:MT,:)=[];
    kwh=kwh';
    kwh=kwh/1000;
    %kwh_num=kwh(bchoice_ij);
    %kwh_denom=kwh(id_ij);
    %clear kwh;
    
    %tmp_file8=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\elec_county_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file8=strcat(path_data{1},'elec_county_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    elec=load(tmp_file8); 
    elec(:,1:7)=[];
    elec(Nsample+1:MT,:)=[];
    elec=elec_scale*elec;

    %elec_J=repmat(elec,1,J);
    %elec_J=elec_J';
    %elec_num=elec_J(bchoice_ij).*kwh_num;
    %elec_denom=elec_J(id_ij).*kwh_denom;
    %clear elec_J;
    
    %tmp_file9=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\demo_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file9=strcat(path_data{1},'demo_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    demo=load(tmp_file9); 

        %hd_id income2 income3 education2 education3 fam_size age political2 political3
    demo(:,1)=[];
    demo(:,1:2)=[];
    demo(Nsample+1:MT,:)=[];
    demo(:,3)=demo(:,3)/10;
    demo(:,4)=demo(:,4)/100;
    demo1_J=repmat(demo(:,1),1,J); %education2
    demo2_J=repmat(demo(:,2),1,J); %education3
    demo3_J=repmat(demo(:,3),1,J); %fam_size
    demo4_J=repmat(demo(:,4),1,J); %age
    demo5_J=repmat(demo(:,5),1,J); %political2
    demo6_J=repmat(demo(:,6),1,J); %political3
   
%Create Interactions Attributes-Demographics     
    %tmp_file10=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\type_id_2_3_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    tmp_file10=strcat(path_data{1},'type_id_2_3_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    type_id2_3=load(tmp_file10); 
    type_id2_3(:,1:7)=[];
    type_id2_3(Nsample+1:MT,:)=[];
   
    type_id2_3_demo1=type_id2_3.*demo1_J;
    type_id2_3_demo1=type_id2_3_demo1';
    type_id2_3_demo1_num=type_id2_3_demo1(bchoice_ij);
    type_id2_3_demo1_denom=type_id2_3_demo1(id_ij);
    clear type_id2_3_demo1;
    
    type_id2_3_demo2=type_id2_3.*demo2_J;
    type_id2_3_demo2=type_id2_3_demo2';
    type_id2_3_demo2_num=type_id2_3_demo2(bchoice_ij);
    type_id2_3_demo2_denom=type_id2_3_demo2(id_ij);
    clear type_id2_3_demo2;
    
    type_id2_3_demo3=type_id2_3.*demo3_J;
    type_id2_3_demo3=type_id2_3_demo3';
    type_id2_3_demo3_num=type_id2_3_demo3(bchoice_ij);
    type_id2_3_demo3_denom=type_id2_3_demo3(id_ij);
    clear type_id2_3_demo3;
    
    type_id2_3_demo4=type_id2_3.*demo4_J;
    type_id2_3_demo4=type_id2_3_demo4';
    type_id2_3_demo4_num=type_id2_3_demo4(bchoice_ij);
    type_id2_3_demo4_denom=type_id2_3_demo4(id_ij);
    clear type_id2_3_demo4;
    
    type_id2_3_demo5=type_id2_3.*demo5_J;
    type_id2_3_demo5=type_id2_3_demo5';
    type_id2_3_demo5_num=type_id2_3_demo5(bchoice_ij);
    type_id2_3_demo5_denom=type_id2_3_demo5(id_ij);
    clear type_id2_3_demo5;
    
    type_id2_3_demo6=type_id2_3.*demo6_J;
    type_id2_3_demo6=type_id2_3_demo6';
    type_id2_3_demo6_num=type_id2_3_demo6(bchoice_ij);
    type_id2_3_demo6_denom=type_id2_3_demo6(id_ij);
    clear type_id2_3_demo6;
   
    clear type_id2_3;
   
    %tmp_file11=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\AV_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    tmp_file11=strcat(path_data{1},'AV_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');      
    AV=load(tmp_file11); 
    AV(:,1:7)=[];
    AV(Nsample+1:MT,:)=[];
    AV=AV/100;
    
    AV_demo1=AV.*demo1_J;
    AV_demo1=AV_demo1';
    AV_demo1_num=AV_demo1(bchoice_ij);
    AV_demo1_denom=AV_demo1(id_ij);
    clear AV_demo1;
    
    AV_demo2=AV.*demo2_J;
    AV_demo2=AV_demo2';
    AV_demo2_num=AV_demo2(bchoice_ij);
    AV_demo2_denom=AV_demo2(id_ij);
    clear AV_demo2;
    
    AV_demo3=AV.*demo3_J;
    AV_demo3=AV_demo3';
    AV_demo3_num=AV_demo3(bchoice_ij);
    AV_demo3_denom=AV_demo3(id_ij);
    clear AV_demo3;
    
    AV_demo4=AV.*demo4_J;
    AV_demo4=AV_demo4';
    AV_demo4_num=AV_demo4(bchoice_ij);
    AV_demo4_denom=AV_demo4(id_ij);
    clear AV_demo4;
    
    AV_demo5=AV.*demo5_J;
    AV_demo5=AV_demo5';
    AV_demo5_num=AV_demo5(bchoice_ij);
    AV_demo5_denom=AV_demo5(id_ij);
    clear AV_demo5;
    
    AV_demo6=AV.*demo6_J;
    AV_demo6=AV_demo6';
    AV_demo6_num=AV_demo6(bchoice_ij);
    AV_demo6_denom=AV_demo6(id_ij);
    clear AV_demo6;
   
    clear AV;
    
    %tmp_file12=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\ice_sc_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    tmp_file12=strcat(path_data{1},'ice_sc_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    ice_sc=load(tmp_file12); 
    ice_sc(:,1:7)=[];
    ice_sc(Nsample+1:MT,:)=[];
    
    ice_sc_demo1=ice_sc.*demo1_J;
    ice_sc_demo1=ice_sc_demo1';
    ice_sc_demo1_num=ice_sc_demo1(bchoice_ij);
    ice_sc_demo1_denom=ice_sc_demo1(id_ij);
    clear ice_sc_demo1;
    
    ice_sc_demo2=ice_sc.*demo2_J;
    ice_sc_demo2=ice_sc_demo2';
    ice_sc_demo2_num=ice_sc_demo2(bchoice_ij);
    ice_sc_demo2_denom=ice_sc_demo2(id_ij);
    clear ice_sc_demo2;
    
    ice_sc_demo3=ice_sc.*demo3_J;
    ice_sc_demo3=ice_sc_demo3';
    ice_sc_demo3_num=ice_sc_demo3(bchoice_ij);
    ice_sc_demo3_denom=ice_sc_demo3(id_ij);
    clear ice_sc_demo3;
    
    ice_sc_demo4=ice_sc.*demo4_J;
    ice_sc_demo4=ice_sc_demo4';
    ice_sc_demo4_num=ice_sc_demo4(bchoice_ij);
    ice_sc_demo4_denom=ice_sc_demo4(id_ij);
    clear ice_sc_demo4;
    
    ice_sc_demo5=ice_sc.*demo5_J;
    ice_sc_demo5=ice_sc_demo5';
    ice_sc_demo5_num=ice_sc_demo5(bchoice_ij);
    ice_sc_demo5_denom=ice_sc_demo5(id_ij);
    clear ice_sc_demo5;
    
    ice_sc_demo6=ice_sc.*demo6_J;
    ice_sc_demo6=ice_sc_demo6';
    ice_sc_demo6_num=ice_sc_demo6(bchoice_ij);
    ice_sc_demo6_denom=ice_sc_demo6(id_ij);
    clear ice_sc_demo6;
   
    clear ice_sc;
      
    prod=repmat([1:J],Nsample,1)';
    prod_denom=prod(id_ij);
    prod_num=prod(bchoice_ij);
    
    
    
    %tmp_file13=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\estarkwh_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file13=strcat(path_data{1},'estarkwh_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    estarkwh=load(tmp_file13); 
    estarkwh(:,1:7)=[];
    estarkwh(Nsample+1:MT,:)=[];
    estarkwh=estarkwh';
    mefkwh=estarkwh/0.8;
    mefkwh=mefkwh/1000;
    mefkwh_num=mefkwh(bchoice_ij);
    mefkwh_denom=mefkwh(id_ij);
    
    %tmp_file14=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\emfac_zip_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    tmp_file14=strcat(path_data{1},'emfac_zip_046_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');    
    em_fac = load(tmp_file14); 
    em_fac(Nsample+1:MT,:)=[];
    d_ext = repmat(d_ext_carbon_price*em_fac',J,1);
    
 clear data; 

  

 	%tmp_file1=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Data\choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    tmp_file1=strcat(path_data{1},'choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
    id=load(tmp_file1); 
	Nmax=size(id,1);
	if subsample == 11000
		bs_max=10*ceil(Nmax/10000);
    elseif subsample == 44000
		bs_max=10*ceil(Nmax/20000); 
    end

c=1
%for j=1:bs_max
for j=1:16
%for j=1:2	
   
	try	
	
	    if subsample == 11000
            index_bs=[(1+(j-1)*(1000)):min(j*(1000),Nmax)];
        elseif subsample == 44000
            index_bs=[(1+(j-1)*(2000)):min(j*(2000),Nmax)];
        end 
        
	    tmp_file1=strcat(path_data{1},'choiceset_identifier_trimester_week_zipcode_2008_2012_v11022017_struct_',num2str(subsample),'_',num2str(inc),'_seed_',num2str(set_seed),'.csv');
	    id=load(tmp_file1);
	    id=id(index_bs,:);
	    id(:,1:7)=[];
	    choicest_size=sum(id,2);
	    %id_discard=find(choicest_size>200 | choicest_size<50);
	    id_discard=find(choicest_size>400);
	    id(id_discard,:)=[];  
	    choicest_size(id_discard)=[]; 
	    id=id';
	    J_bs=size(id,1);
	
		results_file=strcat(path_data{1},'Fox_discrete_mixedFEdemoshrt_WTP_v',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_cty_tau1_grid_v3_','bs_',num2str(j),'.mat');
	    fox=load(results_file);
	    j 
	    eta_fox_v{c}   = fox.param_v(:,1);
	    theta_m_5_v{c} = fox.param_v(:,2);
        tau_fox_v{c}   = fox.param_v(:,3);
	    weight{c}      = fox.beta;
           
        eta_hom(c)   = eta_fox_v{c}'*weight{c};
        theta_hom(c) = theta_m_5_v{c}'*weight{c};
        
        %Adjustment to rescal m:
        theta_m_2_v    = theta_m_5_v{c}*LLC_5_18/LLC_2_18;
        theta_m_12_v   = theta_m_5_v{c}*LLC_5_18/LLC_12_18;
        
        m_2=LLC_2_18/LLC_5_18;
        m_12=LLC_12_18/LLC_5_18;
        
        id_m_2  = (theta_m_5_v{c} > m_2);
        id_m_12 = (theta_m_5_v{c} < m_12);
        id_neg  = (theta_m_5_v{c}<0);
        id_one  = (theta_m_5_v{c}>1);
        
        theta_m_adj = theta_m_5_v{c}*0+1;
        theta_m_adj(id_m_2)  = theta_m_2_v(id_m_2);
        theta_m_adj(id_m_12) = theta_m_12_v(id_m_12);
        theta_m_adj(id_neg)  = theta_m_12_v(id_neg)*0;
        theta_m_adj(id_one)  = theta_m_12_v(id_one)*0+1;
        
        theta_fox_v{c}=theta_m_adj;
        
	   %distribution of parameters used
	    matrix_fox_results{c}=[eta_fox_v{c}, theta_fox_v{c}, weight{c}];
	    
% 	    result_homFE=strcat('\\c3\rdat\SHoude\Research\sears\estar_results\struct_est\hom_FEdemoshrt_wtp_LLC_5_18_',num2str(subsample),'_inc',num2str(inc),'_sd',num2str(set_seed),'_04132018_bs_',num2str(j),'.mat');
% 		eval(['load ' result_homFE]);
% 		result_homFE0=theta_est1;
% 
% 		param_FE{c}=result_homFE0(1:J_bs-1,1);
%  		param_FEdemo{c}=result_homFE0(J_bs+4:J_bs+15)
%  
%  		eta(c) = result_homFE0(J_bs);  
%  		tau(c) = result_homFE0(J_bs+1);
%  	    pi(c) = result_homFE0(J_bs+2);
%  		theta(c) = result_homFE0(J_bs+3);
 		
 		result_mixed=strcat(path_data{1},'mixed_logit_FEdemoshrt_grID_wtp_100iters_',num2str(subsample),'_',num2str(inc),'_',num2str(set_seed),'_matlab_seed_',num2str(matlab_seed),'_nb_draws_',num2str(nb_draws),'bs_',num2str(j),'.mat');
		eval(['load ' result_mixed]);
		result_mixed0=theta_est1;
		param_FE{c}=result_mixed0(1:J_bs-1,1);
		param_FEdemo{c}=result_mixed0(J_bs+4:J_bs+15);
		
  		tau(c)=result_mixed0(J_bs+1);
		pi(c)=result_mixed0(J_bs+2);

		id_eta_neg{c} = find(matrix_fox_results{c}(:,1)<=0);
		eta_restricted{c}=matrix_fox_results{c}(id_eta_neg{c},:);
		eta_mean_experience(c)=eta_restricted{c}(:,1)'*eta_restricted{c}(:,3);
				
        %eta_hom(c)= matrix_fox_results{c}(:,1)'*matrix_fox_results{c}(:,3);
        %theta_hom(c)= matrix_fox_results{c}(:,2)'*matrix_fox_results{c}(:,3);
       
 		id_weight{c} = find(matrix_fox_results{c}(:,3)>=5e-4);
 	    results_k{c}=matrix_fox_results{c}(id_weight{c},:);
 		dim_k{c}=size(results_k{c},1);
 
       eta_theta_k=[];
		 for k=1:dim_k{c}
		   if results_k{c}(k,1)>0
		     eta_experience(c) = eta_mean_experience(c); 
		   else
		     eta_experience(c) = results_k{c}(k,1); 
           end    
		     eta_theta_k{k}=[results_k{c}(k,1),eta_experience(c),results_k{c}(k,2),tau(c),pi(c)];
         end
        weight_c=results_k{c}(:,3);

		%% 
		% Basecase 
		t_pigou=0;
		   
		  [Total_base_SW{c},Total_base_SW_experience{c},Total_base_W_avg{c},Total_base_ext_avg{c},Total_base_gvt_avg{c}] = ...
		  Total_Welfare_clogit_WTP_EEgap_parfor(weight_c,eta_theta_k,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18)
		                                                                         
		Total_base_SW_bs{inc,c}            = Total_base_SW{c};
		Total_base_SW_experience_bs{inc,c} = Total_base_SW_experience{c};
		Total_base_W_avg_bs{inc,c}         = Total_base_W_avg{c};
		Total_base_ext_avg_bs{inc,c}       = Total_base_ext_avg{c};
		Total_base_gvt_avg_bs{inc,c}       = Total_base_gvt_avg{c};  
		
		kwh_denom=kwh(id_ij);  
		cost_base_denom=(191/1000)./kwh_denom;
		
		%%
		% Scenario 1: Optimal Pigouvian Tax 
		
		%parpool(24)
		Total_SW=zeros(1,140);
		Total_SW_experience=zeros(1,140);
		Total_W_avg=zeros(1,140);
		Total_ext_avg=zeros(1,140);
		Total_gvt_avg=zeros(1,140);
		%t_pigou=zeros(1,140);
		u_standard=zeros(1,140);
        
		parfor i=1:140
		 t_pigou=0;
		 u_standard(i)=(100+(i-1)*5)/1000
		 kwh_standard=kwh*0+(100+(i-1)*5)/1000;
		 kwh_standard_denom=kwh_standard(id_ij);
		 
		 cost_standard_denom=(191/1000)./kwh_standard_denom;
		 diff_price_denom=cost_standard_denom-cost_base_denom;
		 price_standard_denom=price_denom+diff_price_denom;
		  
		 [Total_SW(i),Total_SW_experience(i),Total_W_avg(i),Total_ext_avg(i),Total_gvt_avg(i)] = ...
		 Total_Welfare_clogit_WTP_EEgap_parfor(weight_c,eta_theta_k,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_standard_denom, rebate_estar_denom, estar_denom, kwh_standard, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18)
		i
		end
		
		[max_heter, max_id_heter]=max(Total_SW);
		u_standard_star_heter(c)=u_standard(max_id_heter); 
		kwh_standard_star_heter=kwh*0+u_standard_star_heter(c);
		kwh_standard_star_heter_denom=kwh_standard_star_heter(id_ij);
		
		cost_standard_star_heter_denom=(191/1000)./kwh_standard_star_heter_denom;
		diff_price_star_heter_denom=cost_standard_star_heter_denom-cost_base_denom;
		price_standard_star_heter_denom=price_denom+diff_price_star_heter_denom;

		Total_heter_SW_bs{inc,c}            = Total_SW;
		Total_heter_SW_experience_bs{inc,c} = Total_SW_experience;
		Total_heter_W_avg_bs{inc,c}         = Total_W_avg;
		Total_heter_ext_avg_bs{inc,c}       = Total_ext_avg;
		Total_heter_gvt_avg_bs{inc,c}       = Total_gvt_avg;
		
		weight_hom=1
		eta_theta_hom_c{1}=[eta_hom(c), eta_hom(c), theta_hom(c), tau(c), pi(c)];
		
		%%
		% Basecase: Homogeneous M
		t_pigou=0;
		   
		[Total_base_hom_SW{c},Total_base_hom_SW_experience{c},Total_base_hom_W_avg{c},Total_base_hom_ext_avg{c},Total_base_hom_gvt_avg{c}] = ...
		Total_Welfare_clogit_WTP_EEgap_parfor(weight_hom,eta_theta_hom_c,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18)
		                                                                         
		Total_base_hom_SW_bs{inc,c}            = Total_base_hom_SW{c};
		Total_base_hom_SW_experience_bs{inc,c} = Total_base_hom_SW_experience{c};
		Total_base_hom_W_avg_bs{inc,c}         = Total_base_hom_W_avg{c};
		Total_base_hom_ext_avg_bs{inc,c}       = Total_base_hom_ext_avg{c};
		Total_base_hom_gvt_avg_bs{inc,c}       = Total_base_hom_gvt_avg{c};  
		
		%%
		%Optimal policy with homogeneous M 
		
		Total_hom_SW=zeros(1,140);
		Total_hom_SW_experience=zeros(1,140);
		Total_hom_W_avg=zeros(1,140);
		Total_hom_ext_avg=zeros(1,140);
		Total_hom_gvt_avg=zeros(1,140);
		%t_pigou=zeros(1,140);
		u_standard=zeros(1,140);
        
		parfor i=1:140
		 t_pigou=0;
		 u_standard(i)=(100+(i-1)*5)/1000
		 kwh_standard=kwh*0+(100+(i-1)*5)/1000;
		 kwh_standard_denom=kwh_standard(id_ij);
		 
		 cost_standard_denom=(191/1000)./kwh_standard_denom;
		 diff_price_denom=cost_standard_denom-cost_base_denom;
		 price_standard_denom=price_denom+diff_price_denom;
		 
		 [Total_hom_SW(i),Total_hom_SW_experience(i),Total_hom_W_avg(i),Total_hom_ext_avg(i),Total_hom_gvt_avg(i)] = ...
		 Total_Welfare_clogit_WTP_EEgap_parfor(weight_hom,eta_theta_hom_c,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_standard_denom, rebate_estar_denom, estar_denom, kwh_standard, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18)
		 d=i;
		 d
		end

		Total_hom_SW_bs{inc,c}            = Total_hom_SW;
		Total_hom_SW_experience_bs{inc,c} = Total_hom_SW_experience;
		Total_hom_W_avg_bs{inc,c}         = Total_hom_W_avg;
		Total_hom_ext_avg_bs{inc,c}       = Total_hom_ext_avg;
		Total_hom_gvt_avg_bs{inc,c}       = Total_hom_gvt_avg;
		
		 results_file1=strcat(path_results{1},'Optimal_m_adj_pos_Wglow_u_standard_heter_em_fac_EEgap_d0.02_', num2str(subsample),'_inc_',num2str(inc), '_sd', num2str(set_seed),'_bs_',num2str(j),'_FOX_bs.mat');    
		 eval(['save ' results_file1 ' t_pigou Total_SW Total_SW_experience Total_W_avg Total_ext_avg Total_gvt_avg Total_base_SW Total_base_SW_experience Total_base_W_avg Total_base_ext_avg Total_base_gvt_avg Total_base_hom_SW Total_base_hom_SW_experience Total_base_hom_W_avg Total_base_hom_ext_avg Total_base_hom_gvt_avg Total_hom_SW Total_hom_SW_experience Total_hom_W_avg Total_hom_ext_avg Total_hom_gvt_avg']);  
		   
		[max_hom, max_id_hom]=max(Total_hom_SW);
		u_standard_star_hom(c)=u_standard(max_id_hom); 
		kwh_standard_star_hom=kwh*0+u_standard_star_hom(c);
		kwh_standard_star_hom_denom=kwh_standard_star_hom(id_ij);
		
		cost_standard_star_hom_denom=(191/1000)./kwh_standard_star_hom_denom;
		diff_price_star_hom_denom=cost_standard_star_hom_denom-cost_base_denom;
		price_standard_star_hom_denom=price_denom+diff_price_star_hom_denom; 
		
		%Evaluate welfare at different policy	 
		[Total_pigou_heter_SW_heter(c),Total_pigou_heter_SW_experience_heter(c),Total_pigou_heter_W_heter(c),Total_pigou_heter_ext_heter(c),Total_pigou_heter_gvt_heter(c)] = ...
		Total_Welfare_clogit_WTP_EEgap_parfor(weight_c,eta_theta_k,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_standard_star_heter_denom, rebate_estar_denom, estar_denom, kwh_standard_star_heter, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18);
		
		[Total_pigou_heter_SW_hom(c),Total_pigou_heter_SW_hom_experience(c),Total_pigou_heter_W_hom_avg(c),Total_pigou_heter_ext_hom_avg(c),Total_pigou_heter_gvt_avg_hom(c)] = ...
		Total_Welfare_clogit_WTP_EEgap_parfor(weight_hom,eta_theta_hom_c,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_standard_star_heter_denom, rebate_estar_denom, estar_denom, kwh_standard_star_heter, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18);
		
		[Total_pigou_hom_SW_heter(c),Total_pigou_hom_SW_experience_heter(c),Total_pigou_hom_W_heter(c),Total_pigou_hom_ext_heter(c),Total_pigou_hom_gvt_heter(c)] = ...
		Total_Welfare_clogit_WTP_EEgap_parfor(weight_c,eta_theta_k,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_standard_star_hom_denom, rebate_estar_denom, estar_denom, kwh_standard_star_hom, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18);
		
		[Total_pigou_hom_SW_hom(c),Total_pigou_hom_SW_hom_experience(c),Total_pigou_hom_W_hom_avg(c),Total_pigou_hom_ext_hom_avg(c),Total_pigou_hom_gvt_avg_hom(c)] = ...
		Total_Welfare_clogit_WTP_EEgap_parfor(weight_hom,eta_theta_hom_c,param_FE{c},param_FEdemo{c},id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,Nsample,price, price_standard_star_hom_denom, rebate_estar_denom, estar_denom, kwh_standard_star_hom, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18);

	    c=c+1;
	catch
		c=c;	
	end
    
end

Result_base_SW_heter = [zeros(c-1,1),cell2mat(Total_base_SW)',cell2mat(Total_base_SW_experience)',cell2mat(Total_base_W_avg)',cell2mat(Total_base_ext_avg)',cell2mat(Total_base_gvt_avg)'];
Result_base_hom_SW_heter = [zeros(c-1,1),cell2mat(Total_base_hom_SW)',cell2mat(Total_base_hom_SW_experience)',cell2mat(Total_base_hom_W_avg)',cell2mat(Total_base_hom_ext_avg)',cell2mat(Total_base_hom_gvt_avg)'];

Result_t_heter_SW_heter = [u_standard_star_heter' , Total_pigou_heter_SW_heter' ,Total_pigou_heter_SW_experience_heter' ,Total_pigou_heter_W_heter' ,Total_pigou_heter_ext_heter' ,Total_pigou_heter_gvt_heter' ];
Result_t_heter_SW_hom   = [u_standard_star_heter' , Total_pigou_heter_SW_hom' ,Total_pigou_heter_SW_hom_experience' ,Total_pigou_heter_W_hom_avg' ,Total_pigou_heter_ext_hom_avg' ,Total_pigou_heter_gvt_avg_hom' ];
Result_t_hom_SW_heter   = [u_standard_star_hom' , Total_pigou_hom_SW_heter' ,Total_pigou_hom_SW_experience_heter' ,Total_pigou_hom_W_heter' ,Total_pigou_hom_ext_heter' ,Total_pigou_hom_gvt_heter' ];
Result_t_hom_SW_hom     = [u_standard_star_hom' , Total_pigou_hom_SW_hom' ,Total_pigou_hom_SW_hom_experience' ,Total_pigou_hom_W_hom_avg' ,Total_pigou_hom_ext_hom_avg' ,Total_pigou_hom_gvt_avg_hom' ];

Result_base_SW_heter_mean    = full(sum(Result_base_SW_heter,1)/(c-1));
Result_t_heter_SW_heter_mean = full(sum(Result_t_heter_SW_heter,1)/(c-1));
Result_t_heter_SW_hom_mean = full(sum(Result_t_heter_SW_hom,1)/(c-1));
Result_t_hom_SW_heter_mean = full(sum(Result_t_hom_SW_heter,1)/(c-1));
Result_t_hom_SW_hom_mean = full(sum(Result_t_hom_SW_hom,1)/(c-1));

Result_base_SW_heter_se    = full(std(Result_base_SW_heter,1)/((c-1)^.5));
Result_t_heter_SW_heter_se = full(std(Result_t_heter_SW_heter)/((c-1)^.5));
Result_t_heter_SW_hom_se = full(std(Result_t_heter_SW_hom)/((c-1)^.5));
Result_t_hom_SW_heter_se = full(std(Result_t_hom_SW_heter)/((c-1)^.5));
Result_t_hom_SW_hom_se = full(std(Result_t_hom_SW_hom)/((c-1)^.5));

Delta_Result_t_heter_SW_heter_mean = full(sum(Result_t_heter_SW_heter-Result_base_SW_heter,1)/(c-1));
Delta_Result_t_hom_SW_heter_mean = full(sum(Result_t_hom_SW_heter-Result_base_hom_SW_heter,1)/(c-1));
Delta_Result_t_heter_SW_heter_se = full(std(Result_t_heter_SW_heter-Result_base_SW_heter)/((c-1)^.5));
Delta_Result_t_hom_SW_heter_se = full(std(Result_t_hom_SW_heter-Result_base_hom_SW_heter)/((c-1)^.5));

All_results_heter=[Result_base_SW_heter_mean; Result_base_SW_heter_se;...
                   Result_t_heter_SW_heter_mean; Result_t_heter_SW_heter_se;...
                   Result_t_hom_SW_heter_mean; Result_t_hom_SW_heter_se];
               
Delta_results_heter=[ Delta_Result_t_heter_SW_heter_mean; Delta_Result_t_heter_SW_heter_se;...
                   Delta_Result_t_hom_SW_heter_mean; Delta_Result_t_hom_SW_heter_se];  
               
results_file2=strcat(path_results{1},'Tables_Welfare_m_adj_pos_Wglow_Optimal_u_standard_heter_em_fac_EEgap_d0.02_', num2str(subsample),'_inc_',num2str(inc), '_sd', num2str(set_seed), '_FOX_bs.mat');    
eval(['save ' results_file2 ' Delta_results_heter All_results_heter Result_t_heter_SW_heter_mean Result_t_heter_SW_hom_mean Result_t_hom_SW_heter_mean Result_t_hom_SW_hom_mean Result_t_heter_SW_heter_se Result_t_heter_SW_hom_se Result_t_hom_SW_heter_se Result_t_hom_SW_hom_se Total_pigou_heter_SW_heter Total_pigou_heter_SW_experience_heter Total_pigou_heter_W_heter Total_pigou_heter_ext_heter Total_pigou_heter_gvt_heter Total_pigou_heter_SW_hom Total_pigou_heter_SW_hom_experience Total_pigou_heter_W_hom_avg Total_pigou_heter_ext_hom_avg Total_pigou_heter_gvt_avg_hom  Total_pigou_hom_SW_heter Total_pigou_hom_SW_experience_heter Total_pigou_hom_W_heter Total_pigou_hom_ext_heter Total_pigou_hom_gvt_heter Total_pigou_hom_SW_hom Total_pigou_hom_SW_hom_experience Total_pigou_hom_W_hom_avg Total_pigou_hom_ext_hom_avg Total_pigou_hom_gvt_avg_hom']);  
 
results_file3=strcat(path_results{1},'Results_Welfare_m_adj_pos_Wglow_Detailed_u_standard_heter_em_fac_EEgap_d0.02_', num2str(subsample),'_inc_',num2str(inc), '_sd', num2str(set_seed), '_FOX_bs.mat');    
eval(['save ' results_file3 ' inc c Total_base_SW_bs Total_base_SW_experience_bs Total_base_W_avg_bs Total_base_ext_avg_bs Total_base_gvt_avg_bs Total_heter_SW_bs  Total_heter_SW_experience_bs Total_heter_W_avg_bs Total_heter_ext_avg_bs Total_heter_gvt_avg_bs Total_hom_SW_bs Total_hom_SW_experience_bs Total_hom_W_avg_bs Total_hom_ext_avg_bs Total_hom_gvt_avg_bs ']);
		

end 






      












